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Abstract. At low energies, electrons in doped graphene sheets are described by 
a massless Dirac fermion Hamiltonian. In this work we present a semi-analytical 
expression for the dynamical density-density linear-response function of noninteracting 
massless Dirac fermions (the so-called "Lindhard" function) at finite temperature. This 
result is crucial to describe finite-temperature screening of interacting massless Dirac 
fermions within the Random Phase Approximation. In particular, we use it to make 
quantitative predictions for the specific heat and the compressibility of doped graphene 
sheets. We find that, at low temperatures, the specific heat has the usual normal- 
Fermi-liquid linear-in-temperature behavior, with a slope that is solely controlled by 
the renormalized quasiparticle velocity. 
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1. Introduction 

Graphene is a newly realized two-dimensional (2D) electron system that has attracted 
a great deal of interest in the scientific community because of the new physics which it 
exhibits and because of its potential as a new material for electronic technology [U [2]. 
The agent responsible for many of the interesting electronic properties of graphene 
sheets is the non-Bravais honeycomb-lattice arrangement of Carbon atoms, which leads 
to a gapless semiconductor with valence and conduction 7r-bands. States near the 
Fermi energy of a graphene sheet are described by a spin-independent massless Dirac 
Hamiltonian [3] 

Hn=v F (T-p, (1) 

where v p is the Fermi velocity, which is density-independent and roughly three-hundred 
times smaller that the velocity of light in vacuum, and cr = (o~ x ,o~ y ) is a vector 
constructed with two Pauli matrices {cr\i = x,y}, which operate on pseudospin 
(sublattice) degrees of freedom. Note that the eigenstates of TCb have a definite chirality 
rather than a definite pseudospin, i.e. they have a definite projection of the honeycomb- 
sublattice pseudospin onto the momentum p. 

When non-relativistic Coulombic electron-electron interactions are added to the 
kinetic Hamiltonian ([1]), graphene represents a new type of many-electron problem, 
distinct from both an ordinary 2D electron gas (EG) and from quantum electrodynamics. 
The Dirac-like wave equation and the chirality of its eigenstates lead indeed to both 
unusual electron-electron interaction effects [U El El [7J [8] and to unusual response to 
external potentials [91 [10] . 

Within this low energy description, the properties of doped graphene sheets depend 
on the dimensionless coupling constant 

a - = ' (2) 

and on an ultraviolet cut-off A = k c /k F . Here g = g s g v = 4 accounts for spin and valley 
degeneracy, k-p = (Aim/ g) 1 / 2 is the Fermi wave number with n the electron density, and 
k c should be assigned a value corresponding to the wavevector range over which the 
continuum model (CD) describes graphene. For definiteness we take k c to be such that 
Tik 2 = (2ir) 2 /Ao, where Ao = 3^/30^/2 is the area of the unit cell in the honeycomb 
lattice, with ao ~ 1.42 A the Carbon-Carbon distance. With this choice 

The continuum model is useful when k c k F , i.e. when A ^> 1. 

Vafek [11] has recently shown that the specific heat of undoped graphene sheets 
presents an anomalous low-temperature behavior showing a logarithmic suppression 
with respect to its noninteracting counterpart, Cy{T — ► 0)/Cy(0) oc T/ln(T). On 
the other hand, in Refs. [HI [7J we have demonstrated (see also Ref. [8]) that doped 
graphene sheets are normal (pseudochiral) Fermi liquids, with Landau parameters that 
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possess, however, a quite distinct behavior from those of conventional 2D EGs. In this 
work we calculate the Helmholtz free energy JF(T) of doped graphene sheets within 
the Random Phase Approximation (RPA) [T2J, [13] . This allows us to access important 
thermodynamic quantities, such as the compressibility and the specific heat, which can 
be calculated by taking appropriate derivatives of the free energy. We show that, at low 
temperatures, the specific heat of doped graphene, contrary to the one of the undoped 
system [11], has the usual linear- in-temperature behavior, which is solely controlled by 
the renormalized velocity of quasiparticles as in a normal Fermi liquid. 

2. The Helmoltz free energy and the Lindhard response function at finite 
temperature 

The free energy T = T§ + J-'mt is usually decomposed into the sum of a noninteracting 
term JF and an interaction contribution JF int . To evaluate the interaction contribution 
to the Helmholtz free energy we follow a familiar strategy [13] by combining a coupling 
constant integration expression for jF int valid for uniform continuum models (h~ = 1 from 
now on), 

with a fluctuation-dissipation-theorem (FDT) expression [13] for the static structure 
factor, 

S W (g, T) — / duj coth (/3u;/2) W,? (q, u, T) . (5) 



Here v q = 2iTe 2 /(eq) is the 2D Fourier transform of the Coulomb potential and 
{3 = (/cb^) _1 - We anticipate that this version of the FDT (in which the frequency 
integration has to be performed over the real-frequency axis) requires care in handling 
the plasmon contribution to jF int (T) (see discussion below). 

The RPA approximation for jF int then follows from the RPA approximation for 

yWfo u T) — * (0)( ^' T) (6) 

where X {Qi w > T) is the noninteracting density-density response-function, 

d 2 k 1 + ss' cos(6» fcifc+q ) 



X (0) (q,uj,T) = g lim V / 

s,s'=± 



(2< 



2 



x ^F(gfe,s) ~ nF{£k+q,s') 

v + £k,s - £k+ q ,s> + ir} 

Here e^s = sv-pk are the Dirac band energies and n F (e) = {exp[/3(e — /i )] + l}^ 1 is the 
usual Fermi-Dirac distribution function, /io = ^q{T) being the noninteracting chemical 
potential. As usual, this is determined by the normalization condition 

r+oo 

fi (T) = de u(e)n F (e) , (8) 
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where u(e) = ge/(2Trv^) is the noninteracting density of states. For T — > one finds 
[Iq{T) = e F — 7r 2 (T/T F ) 2 /6, where T F = e F /k B is the Fermi temperature. The factor in 
the first line of Eq. ((7j), which depends on the angle #fc,fc+g between k and k + q, describes 
the dependence of Coulomb scattering on the relative chirality ss f of the interacting 
electrons. 

After some straightforward algebraic manipulations we arrive at the following 
expressions for the imaginary [3m x^(g, a;,T)] and the real [9Re X (<7> u i T)) parts 
of the noninteracting density- density response function for uj > 0: 



3m x^Hqj w ) ^) 



9_ 

4n 



^2 \ - uj)q 2 f{v F q, uj) 

a=± { 

x [GP(q,cu,T)-Gt\q,co,T) 
+ Q(uj - v F q)q 2 f(uj, v F q) 



-2 6a 



+ H^(q,u,T) 



(9) 



and 



X<°'(<^,T) = |-W 

a=± ^ 



-2fc B Tln[l + e 



T)] 



+ 0(u; - w F <?) 



x gV(^^g) GW( !)W) T)-^(g lWl !r) 



(a)/ 



+ 6(w F g - u)q 2 f(v F q, u) 



71 



-5 a ,.+H^\q,uj,T) 



Here 



f(x,y) 



'x* 



G ( ^\q,uJ,T) 



du 



u z 



1 



exp 



\v F qu ± uj\ — 2a/i 
2fc B T 



(10) 



'11} 



(12) 



and 



H ( «\q,uj,T) 



du 



IT 



exp 



v F qu ± u\ — 2a/j 
2k B T 



(13) 



These semi-analytical expressions for 3?e x {q^jT) and 3m x (?) <^> ^) constitute 
the first important result of this work. In Fig. [T] we have plotted the static response, 
•Re x^(9j0'^~')> as a function of q/k F for different values of T/T F . The temperature 
dependence of the Lindhard function at finite frequency is instead presented in Fig. [2j 
An illustrative plot of the imaginary part of the inverse RPA dielectric function 



e(q,u) 



T) = 1 — v q x^°\q, oj, T) is reported in Fig. [3j 
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Figure 1. The static response function Sfte x^^g, 0, T) [in units of — v(sf)] as a 
function of qjk-p for three values of < T/Tp < 1. 
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Figure 2. Left panel: the real part of the dynamical response function 3? e x*- ) (g, w, T) 
[in units of —v{e-p)\ as a function of q/k? for w = 2ep and three values of < T/Tp < 1. 
Right panel: same as in the left panel but for the imaginary part. 



2 " 




4 0.0 



-0.5 



•1.0 



■1.5 



-2.0 



1 



4 1 



q/k F 



q/kj 



Figure 3. Left panel: Sm [e~ 1 (q,ui,T)] as a function of q/fcp and uj/ey for a gr = 2 
and T = 0. The red solid line is the plasmon dispersion relation. Right panel: same 
as in the left panel but for T = 0.2 Tp (corresponding roughly to room temperature). 
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The coupling constant integration in Eq. (J4J) can be carried out partly analytically 
due to the simple RPA expression ([6]). We find that the interaction contribution to the 
free energy per particle fmt(T) is given by 



/fa* CO 



N 



d 2 q 



x arctan 



(2tt) 2 
v q Qm X (Qi U i O 



— / coth (/W2) 



- V f 



+ 



d[Ke X {0) (q,uJ,T)] 



^ P i,0 



(14) 



In this equation the first term comes from the smooth electron-hole contribution to 
3m Xp X p 7 while the second term comes from the plasmon contribution; u v \ = oo p i(q, T, A) 
is the plasmon dispersion relation at coupling constant A which can be found numerically 
by solving the equation 1 — XvqSte x (Q> u i O = 0- Note that in a usual 2D EG the 
exchange energy starts to matter little for T of order Tp because all the occupation 
numbers are small and the Pauli exclusion principle matters little. In the graphene case 
however exchange interactions with the negative energy sea remain important as long 
as T is small compared to v-pk c /k-B = TpA. 

The free energy calculated according to Eq. (fT4"]) is divergent since it includes the 
interaction energy of the model's infinite sea of negative energy particles. Following 
Vafek [11], we choose the free energy at T = 0, f(T = 0), as our "reference" free 
energy, and thus introduce the regularized quantity 5f = f(T) — f(T = 0). This 
again can be decomposed into the sum of a noninteracting contribution, 5fo(T — > 
0) = —gep7i 2 (T/T F ) 2 /12, and an interaction-induced contribution 5f iri t(T) = /i n t(0 — 
fmt(T = 0), which can be calculated from Eq. ( TT41) . Numerical results for 5/i nt (T) as a 
function of the reduced temperature T/T F are presented in the left panel of Fig. HI 

The low-temperature behavior of the interaction contribution to the free energy can 
be extracted analytically with some patience. After some lenghty but straightfoward 
algebra we find, to leading order in A, 



SMT - 0) 



rr 



T 



a 



1 - a gr £(a gr )] 



where the function defined as in Eq. (14) of Ref. 

128/(ttV) - 32/(tt 2 x 2 ) + 1/x - h(nx/8), with 



InA + R. T. , (15) 
[6], is given by £(x) = 



h{x) 



2x 3 VT 
1 



arctan 



x" 



vT 



X 



In 



x 



+ Vx T ^T 



for x < 1 



for x > 1 



(16) 



4x 3 a/x 2 — 1 \x — \/x 2 — 1, 
The symbol "R. T." in the l.f.s. of Eq. (|15|) indicates "regular terms", i.e. terms that, 
by definition, are finite in the limit A — > oo. Eq. (TT5T) represents the second important 
result of this work. 
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Figure 4. Left panel: the (regularized) interaction contribution to the free energy 
$fint{T) (in units of the Fermi energy ef) as a function of T/Tp for A = 10 2 . Right 
panel: the specific heat CV(T) (in units of &b) as a function of T/Tp. 



Befor concluding this Section, we remind the reader that in Ref. [6] it has been 
proven that the renormalized RPA quasiparticle velocity v* is given, at weak coupling 
and to leading order in A, by 

^_ = 1+ lyti-Q^Kr)] lnA (17) 

3. The specific heat and the compressibility 

The specific heat can be calculated from the second derivative of the Helmholtz free 
energy, C v = -Td 2 [nSf(T)]/dT 2 [J. Numerical results for Cy(T) as a function 
of temperature are reported in Fig. HI We thus see that Sf int (T — > 0) oc T 2 in 
Eq. (1151) implies a conventional Fermi-liquid behavior with a linear-in-T specific heat. 
Moreover, comparing Eq. (fl"5]) with Eq. (fT7|) we find that the ratio between Cy and its 
noninteracting value Cy ■* is given by 

lim — ^ = — , (18) 

a well-known property of normal Fermi liquids [T21 [T3] . We are thus led to conclude, in 
full agreement with the zero-temperature calculations of the quasiparticle energy and 
lifetime performed in Refs. [6l[7j, that doped graphene sheets are normal Fermi liquids. 
Note that the fact that interactions enhance the quasiparticle velocity [see Eq. ( TT7|) ] 
implies that the specific heat of doped graphene sheets is suppressed with respect to its 
noninteracting value. 

The compressibility can be calculated from the following equation 

11 , d 2 [n6f int (T)] 



n 2 n(T) n 2 Ko(T) On 2 



(19) 



| The second derivative is calculated using the full temperature dependent free-energy of the non- 
interacting system, Sfo(T), and not its analytical expression reported above that is valid only for 
T<T F . 
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Figure 5. The dimensionless ratio k(T)/ko(T) as a function of graphene's coupling 
constant a gr for three values of < T/Tp < 0.2. 

where Kq (T) is the inverse compressibility of the noninteracting system at finite 
temperature. In the low-temperature limit l/[n 2 Ko(T — > 0)] = ne-p/2 + 
gnepTi 2 (T/T F ) 2 /4:8. The dependence of the ratio k(T)/k (T) on a gr and T/T F is shown 
in Fig. [51 

4. Conclusions 

In this work we have presented semi-analytical expressions for the real and the imaginary 
parts of the density-density linear-response function of noninteracting massless Dirac 
fermions at finite temperature. These results are very useful to study finite-temperature 
screening within the Random Phase Approximation. For example they can be used to 
calculate the conductivity at finite temperature within Boltzmann transport theory and 
make quantitative comparisons with recent experimental results in unsuspended [TU [15] 
and suspended graphene sheets [HI [T7] . 

The Lindhard function at finite temperature is also extremely useful to calculate 
finite-temperature equilibrium properties of interacting massless Dirac fermions, such 
as the specific heat and the compressibility. For example, in this work we have been 
able to show that, at low temperatures, the specific heat of interacting massless Dirac 
fermions has the usual normal-Fermi-liquid linear-in-temperature behavior, with a slope 
that is solely controlled by the renormalized quasiparticle velocity. 
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